Novel crystal phase in suspensions of hard ellipsoids 
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We present a computer simulation study on the crystalline phases of hard ellipsoids of revolution. 
For aspect ratios > 3 the previously suggested stretched-fcc phase [D. Frenkel and B. M. Mulder, 
Mol. Phys. 55, 1171 (1985)] is replaced by a novel crystalline phase. Its unit cell contains two 
ellipsoids with unequal orientations. The lattice is simple monoclinic. The angle of inclination of 
the lattice, /3, is a very soft degree of freedom, while the two right angles are stiff. For one particular 
value of /3, the close-packed version of this crystal is a specimen of the family of superdense packings 
recently reported [Donev et al., Phys. Rev. Lett. 92, 255506 (2004)]. These results are relevant for 
studies of nucleation and glassy dynamics of colloidal suspensions of ellipsoids. 
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Classical, hard particles such as non-overlapping 
spheres, rods or ellipsoids are widely used as models for 
granular matter, colloidal and molecular fluids, crystals 
and glasses. Their success — and their appeal — lies in the 
fact that the problem of evaluating a many-body parti- 
tion function is reduced to a slightly simpler, geometrical 
problem, namely the evaluation of entropic contributions 
only. This is an advantage, in particular, for computer 
simulations. Hence one of the first applications of com- 
puter simulations was a study of the liquid-solid phase 
transition in hard spheres [l|. 

In this Letter, we re-examine the high-density phase 
behavior of hard ellipsoids of revolution with short aspect 
ratios. This system has been studied in Monte Carlo sim- 
ulations by Frenkel and Mulder in 1985 [2]. Since then, 
the focus of attention has been on the nematic phase and 
the isotropic- nematic transition JEM Biaxial hard 
ellipsoids have also been studied [y, |7[ . But to our knowl- 
edge, the high-density phases have not been investigated 
further. Knowledge of these phases is relevant for studies 
of elongated colloids in general, and it is crucial for the 
study of nucleation [8[ and glassy dynamics [9[ in hard 
ellipsoids. 

At high densities, Frenkel and Mulder assumed that 
the most stable phase was an orientationally ordered solid 
which can be constructed in the following way: A face- 
centered cubic (fee) system of spheres is stretched by a 
factor x in an arbitrary direction n. This transformation 
results in a crystal structure of ellipsoids of aspect ratio 
x, which are oriented along n. As the transformation is 
linear, the density of closest packing is the same as for 
the closest packing of spheres <p = ir/y/TS ~ 0.7405. Re- 
cently, Donev and co-workers showed that ellipsoids can 
be packed more efficiently if non-lattice periodic packings 
(i.e. packings in which a unit cell contains several ellip- 
soids at different orientations) are taken into account [10[ • 
For unit cells containing two particles, they constructed a 
family of packings which reach a density of <j) = 0.770732 
for aspect ratios larger than \/3- 

We have performed Monte Carlo simulations of hard el- 



lipsoids of revolution with aspect ratios a/b = |, 2, 3, 4, 6 
and found that for large parts of the high density phase 
diagram, the lattice crystal suggested by Frenkel and 
Mulder is unstable with respect to a different crystalline 
phase. Its unit cell is simple monoclinic and contains two 
ellipsoids at an angle with respect to each other. We will 
refer to this new phase as SM2 (simple monoclinic with 
a basis of two ellipsoids). 

Simulations were performed at constant particle num- 
ber A 7 ", pressure P and temperature T. The shape of 
the periodic box was allowed to fluctuate, so that the 
crystal unit cell could find its equilibrium shape. This 
was achieved by implementing the Monte Carlo equiv- 
alent of the simulation method by Parrinello and Rah- 
man [ll| 5 (Hi EH • We constructed the initial solid struc- 
tures by stretching an fee hard sphere crystal along the 
[111] direction by a factor of a/b. Hence the simulations 
were started with crystals identical to the ones studied 
by Frenkel and Mulder. For aspect ratios a/b = |,2,3, 
we simulated eight independent systems, each contain- 
ing N w 1700 particles. Simulations started at P = 
50 k B T/$ab 2 . The pressure was lowered in subsequent 
runs until we observed melting to the nematic phase. 
In each run, equilibration lasted roughly 2 million MC 
sweeps, and was followed by 1.8 to 3.2 million MC sweeps 
for calculating thermodynamic averages. (One MC sweep 
consisted of N attempts to move or rotate a particle and 
one attempt to change the box shape, on average.) For 
aspect ratio a/b = 6, we simulated eight independent 
systems with N = 3072 at a pressure P = 46 k B T/Sab 2 . 
Equilibration and data acquisition summed up to more 
than 3 million MC sweeps. One system with a/b = 4, 
N = 1200 and P = 46 k B T/8ab 2 was first simulated 
sampling only rectangular box shapes for a total of 2.2 
million MC sweeps, and was then simulated including 
non-rectangular box shapes at the same and lower pres- 
sures for several million MC sweeps. 

Particle overlap was checked by a routine [3] which 
uses the Vieillard-Baron [15[ and Perram-Wertheim 
Il7j criteria. In a simulation of a system with a/b = 3 and 
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FIG. 1: Snapshot of the SM2 crystal from different angles, 
with a/b = 3, TV = 1728, and P = 46 k B T/8ab 2 . The online 
version is color coded with respect to orientation. 



TV = 1728, 1 million MC sweeps took about 12 hours of 
CPU time on a 1.8 GHz processor. 

All systems with aspect ratios a/b > 3 left the ini- 
tial fee structure in favor of a simple monoclinic lattice 
with a basis of two ellipsoids (SM2). We will discuss this 
structure for a/b = 3 and return to different aspect ratios 
towards the end of this Letter. 

Figure [T] shows a snapshot of a system in the SM2 
phase with a/b = 3. The color code (online) helps distin- 
guish the two directions of orientation which are present 
in the crystal. In the initial configuration, a stretched- 
fec crystal, all ellipsoids were oriented along the z-axis. 
SM2 is mainly the result of a collective re-orientation. 
The two directions of orientation alternate from layer to 
layer. Layers containing ellipsoids of only one orientation 
are here parallel to the y-axis and form an angle with the 
x-axis. Within layers parallel to the xy-plane, the centers 
of mass of the ellipsoids form a nearly triangular lattice. 
It differs from the initial fee crystal by a slight elongation 
along the x-axis. However, the collective re-orientation of 
the ellipsoids displaced their tips in such a way that they 
now form a rectangular lattice. This can be discerned in 
the top view in Fig. [TJ The tips of the neighboring layers 
interlace. As a result, each ellipsoid now has four nearest 
neighbors above and below, whereas in fee, it had three. 
The total number of nearest neighbors has increased from 
12 to 14, which is indicative of a higher packing efficiency. 

The initial triangular symmetry about the z-axis al- 
lows for two additional, equivalent SM2 configurations, 
which are rotated with respect to the one in Fig. [T] by 
±60°. We observed these possibilities as well; in fact, 
only two out of our eight systems assumed the global 
orientation seen in Fig. [U 

The unit cell of SM2 is shown in Fig. O The open cir- 
cles indicate the centers of the two ellipsoids which form 
the basis. The cell is monoclinic. The yellow (light gray) 




FIG. 2: Unit cell of SM2 with a/b = 3 (color online). The 
open circles indicate the centers of the two ellipsoids which 
form the basis. The cell is monoclinic. f3 exhibits large vari- 
ations. The yellow (light gray) ellipsoid is at the origin, the 
green (dark gray) one is at ^(a+b). The orientations are sym- 
metric about the ac-plane. The parameters in Parts a) and 
b) are average values for TV = 1728 and P = 46 k B T/8ab 2 ; cf. 
Fig.[T] Part c) shows the cell at close packing with (3 = 148.3° , 
where it is an instance of the family of packings introduced 
by Donev et al. |10|]. 



ellipsoid is at the origin, the green (dark gray) one is at 
| (a + b). The orientations are symmetric about the ac- 
plane. The parameters used to produce parts a) and b) 
of Fig. [2] are thermal average values obtained from simu- 
lations with TV = 1728 and P = 46 k B T/8ab 2 ; cf. Fig. ffl 
The cell remained monoclinic even when the pressure was 
lowered down to the melting transition into the nematic 
phase. 

The angle of inclination, /?, relaxes extremely slowly. 
The simulations with TV = 1728 were too slow to equili- 
brate this angle. Therefore we carried out a set of very 
long simulations for a smaller system (TV = 432) with 
initial values of (3 in the range 105° < (3 < 150°. After 
more than 100 million Monte Carlo sweeps, there was 
still no clear evidence for a preferred value of (3. Varia- 
tions of 15° in a single simulation were typical, even at 
P = 46 kBT/8ab 2 . Hence, we expect the shear modulus 
in this degree of freedom to be very small. The other 
two angles, a and 7, were stable at 90°, with fluctua- 
tions of < 1°. The associated shear moduli are much 
larger. The reason for this interesting rheological prop- 
erty is that planes of uniform orientation slide well past 
each other in the c-direction only. In some of the long 
simulations, this led to undulations of the lattice in the 
c-direction to the point of planar defects, which would 
spontaneously heal again. 

To find a lower bound for the maximum density of SM2 
(a/b = 3), we performed simulations sampling only the 
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FIG. 3: Equation of state data for a/b = 3 and TV = 1728. 
The open squares show stretched-fcc data by Frenkel and Mul- 
der [2|], the filled squares the higher-density SM2 phase. Also 
shown are the nematic and fluid branches (triangles and cir- 
cles, respectively). Errors on our data are indicated by the 
size of the symbols. The SM2 curve tends to underestimate 
the density slightly since (3 was not equilibrated entirely. 



unit cell parameters and particle orientations, and im- 
posing all symmetries of SM2. The initial parameters 
were average values obtained from the simulations with 
TV = 1728 and P = 46 k B T/8ab 2 . In the process of maxi- 
mizing the density, (3 increased from 105° to about 150°, 
and the common tilt of the ellipsoids with respect to the 
froplane disappeared. We then imposed that (a + c) be 
perpendicular to c (see Fig. [2] Part c)) which is equiva- 
lent to (3 « 148°. Under this condition SM2 becomes an 
instance of the family of packings introduced by Donev 
et al. This simulation achieved the highest packing frac- 
tion, namely <p « 0.770732 (the value reported by Donev 
et al.). 

But already at (3 « 105° we found a jamming density 
of 99.663% of the maximum. Simulations at intermediate 
values indicate a smooth approach towards the maximum 
density as /3 increases. Thus, the close-packing density 
varies very weakly for 105° < (3 < 148°. While this range 
is traversed, ellipsoids of one orientation move past neigh- 
bors of the other orientation by almost half their length. 
This can be seen in Parts b) and c) of Fig. [2l For reasons 
of symmetry, this translation may even continue by the 
same amount while the density remains above 99.663% of 
the maximum. These observations are in accord with the 
fact that (3 is soft at finite pressures; evidently, the free 
volume distribution possesses a similarly slight variation 
with (3. 

Figure [3] shows equation of state data of SM2 from our 
simulations with TV = 1728 particles. The density of SM2 
is higher than that of stretched fee for all pressures. Five 
of eight systems underwent the transition to SM2 already 
at the highest simulated pressure P = 48 k B T/8ab 2 , the 
remaining three at P = 46 k B T/8ab 2 . Note also that in 



all our runs, SM2 melted to the nematic phase without 
re-visiting the stretched fee phase from which it devel- 
oped; evidently, SM2 does not only pack more efficiently 
than stretched fee, it also provides for a better distribu- 
tion of free volume at all densities until the transition to 
the nematic. Hence it is more stable than stretched fee. 
We also show the nematic branch from an (TV, P, T) com- 
pression (i.e. the pressure was raised between successive 
simulations) with TV « 780 particles and up to 6 million 
MC sweeps per run. Even at strong over-compression, 
no spontaneous crystallization occurred. This indicates 
that the nucleation barrier to the SM2 phase is very high. 
Also shown is the isotropic fluid branch as obtained from 
(TV, P, T) compression and expansion runs with TV « 780. 

All eight simulations at a/b = 6 and P = 46 k B T/8ab 2 
formed SM2 as well, although four of them retained a 
planar defect. Different regions in the periodic box were 
able to develop different global orientations of SM2 as the 
systems were larger (TV = 3072) than those with a/b = 3 
(TV = 1728). We also simulated a system with a/b = 4, 
TV = 1200 and P = 46 k B T/8ab 2 ; it formed SM2 as 
well. It also developed a planar defect, this time owing 
to a geometrical mismatch between the simulation box 
and the SM2 unit cell. Note that for a/b = 3 and 4, SM2 
formed even in simulations sampling only rectangular box 
shapes. It is therefore more stable than stretched fee even 
when it cannot assume its equilibrium shape. 

By contrast, ellipsoids with a/b = 2 and the oblate 
a/b — \ showed no tendency to leave the initial stretched- 
fcc structure. We studied each of these systems with eight 
independent simulations. In none of them two preferred 
directions of ellipsoid orientation developed. All of them 
melted to the nematic phase on expansion, directly from 
stretched fee. But note that the apparent stability of 
fee in our simulations may well be due to a free energy 
barrier, rather than indicating genuine stability. 

In Fig. [H we show a phase diagram of hard ellipsoids 
of revolution. It includes part of the results of Frenkel 
and Mulder, and their suggested phase boundaries and 
coexistence regions. We have inserted our state points 
and a vertical dashed line to delimit the region in which 
we found SM2; but our data is insufficient to locate a 
phase boundary. 

In the high density phase diagram of hard ellipsoids of 
revolution we have found a crystal which is more stable 
than the stretched fee structure proposed by Frenkel and 
Mulder [2[. The new phase, SM2, has a simple monoclinic 
unit cell containing a basis of two ellipsoids. The angle 
of inclination (3 is a very soft degree of freedom, whereas 
the other angles are not. At one value of (3 (« 148.3° for 
a/b = 3), close-packed SM2 is an instance of the family of 
packings introduced by Donev et al. [10[. As for thermo- 
dynamic stability, our results unequivocally remove the 
stretched fee structure for aspect ratio a/b = 3 from the 
phase diagram of hard, uniaxial ellipsoids. Our state 
points for a/b = 4 and 6 suggest that this holds for 
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FIG. 4: Phase diagram of hard, uniaxial ellipsoids, show- 
ing the results of Frenkel and Mulder (open symbols) [2], 
and their suggested phase boundaries and coexistence regions. 
The data points at a/b = 1 are taken from [l8[. We have in- 
serted our state points (filled symbols) and a vertical dashed 
line to delimit the region in which we found SM2; but our 
data is insufficient to locate a phase boundary. 



the entire range of 3 < a/b < 6, and possibly beyond. 
However, this does not prove that SM2 is the ground 
state. A procedure has been developed for making al- 
most monodisperse ellipsoids pjj, [20], which are of col- 
loidal size. Their behavior at water-air interfaces has 
been studied [U [22[ ; also 3D structural properties of a 
sedimentation of these particles have been successfully 
characterized [23]. It would be interesting to perform 
experiments probing colloidal crystals of ellipsoids. 
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